%% Meager (2019) benchmark analysis
% Replicates Figure 8 left and center panels, Figure 9, Figure A-10

%% Benchmark grouping
rng(11)

% Keep the indices that has data
profit_indices = [1, 2, 3, 4, 5, 6, 7];
revenue_indices = [1, 2, 3, 4, 5, 6, 7];
expend_indices = [1, 2, 3, 4, 5, 6, 7];
consump_indices = [1, 2, 3, 4, 5];
durables_indices = [2, 3, 4, 5];
tempt_indices = [1, 2, 3, 4, 5];

% Randomly divide studies into groups
group_profit = random_pair(profit_indices);
group_revenue = random_pair(revenue_indices);
group_expend = random_pair(expend_indices);
group_consump = random_pair(consump_indices);
group_durables = random_pair(durables_indices);
group_tempt = random_pair(tempt_indices);

%% Within-group assignment
seed = 11;

% Read the "Estimates" and "se" tables
file_path = '../Data/Meager_Data.xlsx';
estimates = readmatrix(file_path, 'Sheet', 'N=15', 'Range', 'B3:D17');
se = readmatrix(file_path, 'Sheet', 'N=15', 'Range', 'H3:J17');

% Randomly assign baseline and validation studies
[hat_theta, hat_vartheta, rep_base_se2, rep_val_se2] = random_assignment(estimates, se, seed);

%% Case 1: nu=0
rng(11)
N = 15;
sim_cov = 5000;
sim_S = 5000;
nu = 0;

% Compute posterior means and variances for tau
[bar_tau_1, bar_V_tau_1] = compute_posterior_moments(hat_theta, rep_base_se2, nu);

% Check whether the validation study estimates are located in the predicted interval
cov_stat_sim = simulate_cov_stat_distribution(N, sim_cov);
[predict_se_1, empirical_cov_1, cov_p_value_1] = check_CI(hat_vartheta, rep_val_se2, bar_tau_1, bar_V_tau_1, nu, cov_stat_sim);
[interval_1, order] = plot_intervals(hat_vartheta, bar_tau_1, predict_se_1);
filename = sprintf("../Results/nu_%.f_Interval_N_15.png", nu);
saveas(interval_1, filename)

% Generate the S statistics and p-value
S_stat_sim = simulate_S_stat_distribution(N, sim_S);
yticks = 0:0.05:0.35;
[S_stat_1, S_p_value_1, PIT_plot_1] = PITs(hat_vartheta, bar_tau_1, predict_se_1, S_stat_sim, 5, yticks);
filename = sprintf("../Results/nu_%.f_PIT_N_15.png", nu);
saveas(PIT_plot_1, filename);

%% Case 2: EB nu
% Estimate nu
tol = 10^-6;
EB_est_nu = MLE_nu(hat_theta, rep_base_se2, tol);

% Compute posterior means and variances for tau
[bar_tau_2, bar_V_tau_2] = compute_posterior_moments(hat_theta, rep_base_se2, EB_est_nu);

% Check whether the validation study estimates are located in the predicted interval
[predict_se_2, empirical_cov_2, cov_p_value_2] = check_CI(hat_vartheta, rep_val_se2, bar_tau_2, bar_V_tau_2, EB_est_nu, cov_stat_sim);
interval_2 = plot_intervals_fixord(hat_vartheta, bar_tau_2, predict_se_2, order);
filename = sprintf("../Results/EB_Interval_N_15.png");
saveas(interval_2, filename)

% Generate the S statistics and p-value
yticks = 0:0.05:0.35;
[S_stat_2, S_p_value_2, PIT_plot_2] = PITs(hat_vartheta, bar_tau_2, predict_se_2, S_stat_sim, 5, yticks);
filename = sprintf("../Results/EB_PIT_N_15.png");
saveas(PIT_plot_2, filename);

% Plot the ratio of nu to se
nu_se_fig = plot_nu_se_ratio(EB_est_nu, rep_base_se2, rep_val_se2);
saveas(nu_se_fig, "../Results/nu_se_ratio_N_15.png")

%% Report the results
empirical_cov = [empirical_cov_1, empirical_cov_2]
cov_p_value = [cov_p_value_1, cov_p_value_2]
S_stat = [S_stat_1, S_stat_2]
S_p_value = [S_p_value_1, S_p_value_2]
